% -----------------------------------------------------------------------
% This file is part of the ASTRA Toolbox
% 
% Copyright: 2010-2018, imec Vision Lab, University of Antwerp
%            2014-2018, CWI, Amsterdam
% License: Open Source under GPLv3
% Contact: astra@astra-toolbox.com
% Website: http://www.astra-toolbox.com/
% -----------------------------------------------------------------------
addpath(genpath('E:\downloads\astra-2.1.0')) 
addpath(genpath('E:\downloads\NIfTI_20140122')) 

clear all
close all
LD = 430.08;   %1024*0.42=430.08
imN = 1024;
imNZ = 1024;
det_width = 0.42; 
Col_det_count = 1024;     %探测器列数  宽
Row_det_count = 1024;   %探测器行数    高
source_origin = 450;    %SOD射线源-物体距离
origin_det =705-source_origin ;   %DOD物体-探测器距离
N = 360;  %可能需要调整

%% 

FOV = LD * source_origin /( 2* sqrt(LD^2/4 + (source_origin+origin_det)^2)); 
min_x = -FOV;
max_x = FOV;
min_y = -FOV;
max_y = FOV;
min_z = -FOV; 
max_z = FOV;
vol_geom = astra_create_vol_geom(imN, imN, imNZ, min_x, max_x, min_y, max_y, min_z, max_z);

%% 
angles=linspace2(0,2*pi,N);%360度顺时针是2，360度逆时针是-2，180度顺时针是1
proj_geom = astra_create_proj_geom('cone', det_width,det_width,Row_det_count, Col_det_count, angles, source_origin, origin_det);

%% 加载投影数据
addpath('E:\data\1101-4-360tif\')
finf = dir('E:\data\1101-4-360tif\*.tif');   % 数据格式tif
Im= imread (finf(1).name);
%fid=fopen(finf(1).name,'r'); 
%A=fread(fid,[1024,1024],'uint16');
%Im=A'; 
%Im=A
%fclose(fid);

%I1=rgb2gray(Im);
I1=im2gray(Im);
%I1=Im
UnSam=1;
Row = size(I1,1)/UnSam;     %行
Col = size(I1,2)/UnSam;       %列 
Dn = length(finf);    %扫描分度数
projections = zeros(Row,Col,Dn);
for j=1:Dn
      image_name=finf(j).name;             
      Im1=imread(image_name);

      %fid=fopen(image_name,'r');
      %A=fread(fid,[1024,1024],'uint16');
      %Im1=A';
      %Im1=A
      %imtool(A,[]);
      %fclose(fid);



      %I = rgb2gray(Im1);

      I = im2gray(Im1);
      %I=Im1
      I=single(I);
      Ic = imresize(I,[Row Col]);
      projections(:,:,j)=Ic;
end
tic 
projections=log(projections);
projections=max(projections(:))-projections;
proj_data = projections;
proj_data = permute(proj_data, [2 3 1]);
size (proj_data) 

proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, proj_data);

% Create a data object for the reconstruction
rec_id = astra_mex_data3d('create', '-vol', vol_geom); 

% Set up the parameters for a reconstruction algorithm using the GPU
cfg = astra_struct('FDK_CUDA');%SIRT3D_CUDA  FDK_CUDA
cfg.ReconstructionDataId = rec_id;
cfg.ProjectionDataId = proj_id;
% Create the algorithm object from the configuration structure
alg_id = astra_mex_algorithm('create', cfg);
% Run 150 iterations of the algorithm 
astra_mex_algorithm('iterate', alg_id, 150);
% Get the result 
rec = astra_mex_data3d('get', rec_id); 
size(rec)
figure, imshow(squeeze(rec(560,:,:)),[]);   


nii = make_nii(rec); % 形成nii格式的数据
save_nii(nii,strcat('E:\data/','case_6','.nii')); 
% 
% for i=1:size(rec,3)
%     img=rec(:,:,i);
%     metadata=dicominfo('CT-MONO2-16-ankle.dcm');
%     metadata.Width=1024;
%     metadata.Height=1024;
%     ffl=[num2str(i,'%04d'),'.dcm'];
%     dicomwrite(img,ffl,metadata);
% end    



%rec(256,:,:)
% Clean up. Note that GPU memory is tied up in the algorithm object,
% and main RAM in the data objects.
astra_mex_algorithm('delete', alg_id);
astra_mex_data3d('delete', rec_id);
astra_mex_data3d('delete', proj_id);
toc
